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Abstract 

We consider the problem of accurately computing low Mach number flows, with 
the specific intent of studying the interaction of sound waves with incompressible 
flow structures, such as concentrations of vorticity. This is a multiple time (and/or 
space) scales problem, leading to various difficulties in the design of numerical 
methods. In this paper, we concentrate on one of these difficidties - the develop- 
ment of boundary conditions at artificial boundaries which allow sound waves and 
vortices to radiate to the far field. Nonlinear model equations are derived based 
on assumptions about the scaling of the variables. We then linearize these about a 
uniform flow and systematically derive exact boundary conditions using transform 
methods. Finally, we compute useful approximations to the exact conditions which 
are valid for small Mach number and small viscosity. 


1 Introduction 


The generation and interaction of sound waves with complex fluid flows is of great interest 
both from the point of view of fundamental fluid mechanics and applied mathematics 
and from the point of view of practical applications. Many important examples of these 
phenomena occur at low Mach number. A short list includes sound generation and 
propagation in water and the interaction of sound waves with laminar flames. Some 
advantages of studying low Mach number flows are the absence of shocks and the clear 
separation of ‘incompressible’ flow features, namely vortex dynamics, and the sound field. 
A disadvantage is the multiple time scales, requiring accurate computation for very long 
times as measured on the ‘fast’ scale. As a result, the numerical analysis of low Mach 
number flow is somewhat undeveloped in comparison with analytical theories. 

Among the particular difficulties in the construction of accurate and efficient numeri- 
cal methods are the time stepping scheme and the choice of radiation boundary conditions 
at artificial boundaries. The goal in choosing the time stepping procedure is to somehow 
exploit the simplicity of the equations governing the fast dynamics (essentially the wave 
equation) to allow large time steps. The problem with boundary conditions, which is the 
main topic of this work, is to find computationally usable procedures (for example local 
in time conditions) which are accurate over long times. 

Essentially all problems in aero or hydroacoustics involve the radiation of energy to the 
far field in the form of sound waves and vortices. Any computational study then requires 
the introduction of artificial boundaries where accurate boundary conditions must be 
imposed. There is a vast literature concerning the construction of such conditions. (See, 
e.g. [4] and the references contained therein.) This literature is essentially divided into 
separate streams treating either hyperbohc problems or dissipative problems such as 
advection-diffusion equations. Our equations involve the coupling of a hyperbolic system 
governing the sound waves and an advection-diffusion equation for the vorticity, again 
acting on different time scales. 

The remainder of the paper is organized as follows. In Section 2 we present the scal- 
ings we assume hold for the isentropic Navier-Stokes system and use them to derive a 
somewhat simplified set of nonlinear model equations. These are then linearized about 
a uniform flow for the purpose of systematically deriving exact and approximate (in 
the small Mach number - large Reynolds number limit) boundary conditions. The con- 
struction is carried out in Section 3. The result is a reasonably simple set of boundary 
conditions which combine ‘standard’ boundary operators of hyperbolic and advection- 
diffusion type. 
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2 Scalings and Model Equations 

Consider the Navier-Stokes equations for conservation of mass ajid momentum in a fluid 
at constant entropy in two space dimensions: 


Pt -1- upx + VPy -1- pUx + pVy = 0, 

(2.1) 

1 p /4 1 \ 

Ut + UUx -1- VUy -1- -Px = - + “vv + 3^a:v j , 

(2.2) 

1 p /4 1 \ 

Vt -1- UVx + VVy -1- -Py = - + “^XX + ■^'fj-xy J ■ 

(2.3) 


These must be supplemented by an equation of state relating pressure and density. For 
simplicity in the present work we use the equation of state for a 7-law gas: 


p = Kp\ 


(2.4) 


though eventually other equations of state, for example ones appropriate for liquids, will 
be included in the model. 

To put the equations in nondimensional form we introduce a characteristic length 
scale, L, fluid velocity, S, and fluid pressure, P. From the latter we deduce a characteristic 
density, D = and sound speed, — 'yPjD. We then have two important 

dimensionless parameters, the Mach number which we take to be small and the Reynolds 
number which we take to be large: 

M = ^<1, Pe = — >1. (2.5) 

C fi 


There are two natural time scales, based respectively on the fluid velocity and sound 
speed: 

Tslow = ■^) ^fast ^slow (2-6) 

As we are interested in flows where ‘fast’ sound waves are present, we choose the latter 
as our characteristic time scale. We use the same letters for the dimenionless variables 
as for the dimensioned variables above, except for the fast time vaxiable which we call 
T to maintain notational consistency with our study of the one- dimensional problem [7]. 
(We reserve t for the slow time, t = Afr.) We then have: 

Pr + MpUx -f MpVy Mupx + Mvpy = 0, (2.7) 


Ur + Muux + MvUy + 


1 

')M p 





+ 




( 2 . 8 ) 
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Vr + MuVx + MvVy + 


1 

P = 


M 

pRe 



Vxx ”t” 




(2.9) 

( 2 . 10 ) 


A glance at the momentum equations reveals potentially large 0{M ~^ ) terms involv- 
ing the pressure gradient. If indeed these terms were of that order we would expect the 
velocities to become large and the local Mach numbers to become 0(1). Therefore, if 
the Mach number is to be low throughout the flow, the pressure gradient must be 0(M). 
This leads us to introduce a new variable, q, which contains the pressure variations. That 
is; 

p = 1 -|- ‘jMq, p = (1 -t- 'yMqY^'* = 1-1- Mq O(Af^). (2.11) 

Substituting this into the dimensionless system and discarding terms 0{M^) we finally 
obtain our nonlinear model system: 


qr + {l + Mq)ux + (1 + Mq)vy -|- Muq^ + Mvqy = 0, (2-12) 


Ut -F (1 - Mq)qx + Muux + MvUy = — 

M 

Ut + (1 — Mq)qy -t- MuVx + MvVy = — 



d" '^yy 

+ Vxx 


+ ^Vxy) 


J 


(2.13) 

(2.14) 


Alternatively, we could have fixed the time scale, T, and then chosen from two natural 
spatial scales, CT and ST. This suggests the possibility of multiple spatial scales present 
in the solution, which is the case for many important aeroacoustic phenomena. To justify 
the approximations, we must assume that L is chosen so that derivatives are 0(1). For 
certain problems, for example the aeolian tones produced by flow past a cylinder, this 
implies that the sound waves will be slowly varying in space. We have not tried to make 
use of this in our derivation of boundary conditions. 

It is of interest to compare this model system to the equations considered by other 
authors. Both Klainerman and Majda [9] and Kreiss, Lorenz and Naughton [11] have 
studied the incompressible limit, M — > 0. Then it is natural to take the slow time scale, 
where dj dr is replaced by M-dfdt, and to suppose that pressure variations scale like M^. 
Then one obtains the incompressible Navier-Stokes equations by setting M = 0. This 
is a singular perturbation problem, as the incompressible equations require fewer initial 
conditions and fewer boundary conditions at inflow than the compressible equations. 
Hence, there is a possibility of boundary layers for small M, as analyzed in [11]. It is 
important that the conditions we develop do not generate such layers. Requirements on 
the initial data so that these ‘nearly incompressible’ scalings are maintained are studied 
both numerically and theoretically in [3]. 
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Setting M = 0 in our case leads instead to a linear symmetric hyperbolic system 
governing the sound waves: 


9 '' 

u 

/ r 


+ 


0 

1 

\o 


1 0 
0 0 
0 0 


/ 9 \ 
u + 

V ^ 


0 

0 

\1 


0 

0 

0 


1 

0 

0 



(2.15) 


Our boundary conditions must also be accurate for this system, which is equivalent to 
the wave equation for q and the dilatation, Ux + Vj,, coupled with the condition that 
vorticity be constant. (For finite M this means that vorticity is slowly varying.) In the 
next section we shall see that our boundary condition construction involves a standard 
(if unsatisfactorally solved) problem in the design of boundary conditions for the wave 
equation. 

Finally, in order to carry out an analytic study of boundary conditions using transform 
methods, we linearize our model problem about a uniform flow, 17, V . Again keeping the 
same letters for the linearized variables we have: 


qr + Ux + Vy MU qx + MV qy = 0, 

(2.16) 

M 

Ut + qx + MUux + MVuy = — 

^ 2 '^xx “1" '^yy “1" 2 '^xy^ j 

(2.17) 

M 

+ 9y + MUvx -1- MVVy = — 1 


(2.18) 


3 Boundary Conditions for the Model Problem 

To explain our principle of deriving exact boundary conditions, we consider a general 
constant coefficient Cauchy problem of the form 

Wt = ^ 

w{x,y^0) = f{x,y). (3.2) 

It is assumed that f E and that (3.1)-(3.2) is well-posed in L^. (See [10] for a definition 
and discussion of this concept.) Also, we assume the initial data have compact support 
in X. More precisely, we assume that f(x, y) is only nonzero for 

— L-j-6<x<L — S where L > S > 0. (3-3) 

We want to replace the Cauchy problem by an initial-boundary value problem (IB VP) 
with boundary conditions at x = iX so that the solution of the IBVP agrees with the 
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solution of the Cauchy problem restricted io -L < x < L. Applying Fourier transfor- 
mation in y and Laplace transformation in t to (3.1)-(3.2) we arrive at 

sw{x,k,s) — P{—,ik)w{x,k,s) = f{x,k), x G IR. (3-4) 

For any fixed k, s, equation (3.4) is an ODE system in x, which is generally of mixed order. 
We want to derive boundary conditions at x = iL which determine the L^-solution of 
(3.4). To this end, we assume that (3.4) can be written as a first-order system 

-^W - M{k, s)W = F(x, k), X G IR. (3.5) 

ax 

Well-posedness of (3.1)-(3.2) in implies that the symbol P{iki,ik 2 ) of the partial 
differential operator P obeys an estimate 

||gP(ifci.ifc2)t|| < fci,A:2G]R, t>0. (3.6) 

One can then prove that the system matrix M = M(fc, s) in (3.5) has no purely imaginary 
eigenvalues for Re s > ct. Consequently, for Re s ^ a, the system can be transformed 
to a block form in which the exponentially growing modes (in x) are separated from 
the exponentially decaying ones. This is the main idea to obtain the exact boundary 
conditions. We assume that M has dimension D x D with eigenvalues in the left 
half-plane and D+ = D — D- eigenvalues in the right half-plane. We further assume, for 
simplicity, that M can be diagonalized. A nonzero row vector ip of dimension D is called 
a left eigenvector oi M li ip M = Xip. 

Theorem. Let Res > a and assume that M{k,s) has a complete set of eigenvectors. 
Suppose that 


ip^ = ip^{k,s), j = (3.7) 

are linearly independent left eigenvectors of M{k,s) corresponding to eigenvalues with 
negative real parts. Similarly, suppose that 

ip^ = ip^{k,s), j = D- + 1, . . . ,D, (3.8) 

are linearly independent left eigenvectors of M{k,s) corresponding to eigenvalues with 
positive real parts. The boundary conditions 

iP^{k,s)W = 0 atx = -L, i = l,...,D_, (3.9) 

iP^{k,s)W=^0 atx = L, j = D. + (3.10) 
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determine a unique solution of the system (3.5), namely the unique solution in Zr^(lR). 
Inversion of the Fourier-Laplace transform leads to the solution of the Cauchy problem. 

We now apply the general theory to the model system (2.16)-(2.18) where we cissume 
U > 0. (Hence, x = —L corresponds to an inflow boundary and x = L to em outflow 
boundary of the underlying uniform flow.) Since the exact boundary conditions are 
independent of the initial data, the data are ignored in the following. Fourier-Laplace 
transformation leads to 


aq + Ux + ikv + MU qx + ikMV g = 0 , 

( 3 . 11 ) 

/4 , , 2 - ik ^ \ 

au + + MUux + ikMVu = Mvy-^Uxx — « "U + —Vxj, 

( 3 . 12 ) 

av + ikq + MUvx + ikMVv = Mv( — -k^v + Vxx + ^Ux)- 

( 3 . 13 ) 

Here v = If Re, and a = Ms is the dual to the stetched time variable, 
determine the dispersion relation, we make the usual ansatz 

r = t/M. To 


( 3 . 14 ) 

and obtain the condition 

det A(A) = 0 , 

( 3 . 15 ) 

the coefficients of the 3x3 matrix A(A) being at most quadratic in A. The polynomial 
det -A(A) has degree 5 and factors into the product of the quadratic 

Q^{\) = a + ikMV + Muk^ + MUX - MvX^, 

( 3 . 16 ) 

and the cubic 


g 3 (A) = A^ - - (a + ikMV + + MUX - ^A"*) (a + ikMV + 

' o o 

mux). ( 3 . 17 ) 

(This factorization is not accidental. In fact, if one writes (2.16)-(2.18) in terms of 
uj = Vx — Uyj 6 = Ux + Vyj and g, then one obtains a second-order equation for u; which is 
decoupled from a system for (q,S), The quadratic equation ^ 2 (-^) = 0 is the dispersion 
relation for u whereas ( 53 (A) = 0 is the dispersion relation for the (g,^) system.) 

Henceforth we assume 

£7 = 0 ( 1 ), k = 0 ( 1 ). 

( 3 . 18 ) 


Then the solutions ofQ 2 (^) = 0 are 
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with 

U f , 2 vcr .\ 

(3.19) 


a = G + ikMV. 

(3.20) 

Setting M = 0 in Q 3 , one 

obtains the roots 



As, 4 « ±-s/ct2 ^ 

(3.21) 

The third root of Qz is 

^ 3 a ZU 

~ AM^uU MU^ Av' 

(3.22) 


To determine the exact conditions, we need to rewrite (3.11)-(3.13) as a first-order 
system for 


W = {q,u,Ux,v,Vx)^ (3.23) 

of the form 

W, = M{k,s)W, (3.24) 

and then to determine the left eigenvectors of M. The eigenvalues of M are precisely the 
roots Ai, . . . , A 5 . They satisfy the sign conditions: 

ReAj>0 for 7 = 1,3 (3.25) 


ReAj<0 for j = 2,4,5 (3.26) 

To leading order (recall the assumption (3.18)) one obtains for the corresponding left 
eigenvectors of M: 

fs ^ik,0,0,a, MU^, (3.27) 

^ (ik,ikMU,0,cr,0), (3.28) 

-03 « ( - + fc2, a, 0, -ikMU, 0), (3.29) 

0“ « (Vcr 2 + 0 , -ikMU, 0), (3.30) 

0 ® « (a, -MU a, l,ik, 0 ) . (3.31) 
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The corresponding boundary conditions for j = 1,2,5 can be transformed easily to 
physical space. One obtains (to the order given): 

At X = L (from j = 1): 

qy + Vr + MVvy + MUv^ = 0. (3.32) 

At X = —L (from j = 2): 

qy + MUuy + Vr + MVvy = 0. (3.33) 

At X = —L (from j = 5): 

qr + MVqy — MUur + + Uj, = 0 (3.34) 

To translate the conditions for j = 3,4 into physical space, locally in time, we must 
approximate the root -y/cr^ + by a rational function in a. This difficulty, which is 
not unexpected, is precisely the one encountered in developing bovmdary conditions for 
the wave equation. Most approximations proposed in the literature are designed to be 
accurate in the (k/a) — > 0 limit. However, we also want accuracy over long times, 
suggesting the use of an approximation accurate as a = Ms —* 0. Recently, a number 
of papers have appeared dealing with the long time behavior of approximate boundary 
conditons for the wave equation. Barry, Bielak, and MacCamy [1] introduce a useful 
notion of dissipativity, where the large a approximation is modified to avoid exponential 
error growth. Engquist and Halpern [2] propose conditions which are exact in both the 
(7 — > oo and a — ^ 0 limits, and prove that solutions satisfying these rapidly converge to 
the correct steady state. However, it has also been shown [6] that, for certain exterior 
problems in two space dimensions, good long time accuracy is impossible to attain with 
standard boundary conditions. 

A simple approximation, appearing already in [2], is 

+ jfc2 ^a+\k\. (3.35) 

If (J^g){k) denotes the Fourier transform of g{y), we define 

Hg^J^-\\k\{J^g){k)). (3.36) 

Then the conditions for j = 3, 4 translate to: 

At X = L (from j = 3, y/a'^ k^ « cr + |fc|): 

-qr- Hq + ur + MVuy - MUvy = 0. (3.37) 

At X = —L (from j = 4, \/cr 2 + P ~ or -|- |A:|): 

qr + Hq + Ur + MVuy — MUvy = 0. (3.38) 
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Note that the conditions are still nonlocal in y since H is not a local operator. Often, this 
does not lead to difficulties. For example, if the given problem is periodic in y and one 
uses a discrete Fourier method in j/-direction, then the operator H is easily discretized. 
Instead of the approximation (3.35) one can also try 

+ „ > 0. (3.39) 

cr + a 

which goes over into (3.35) for a -> oo. The approximate eigenvectors are multiplied 
through by (7 + o, and one obtains in physical space: 

ki X = L (from j = 3, (3.39)): 

- aHq + (^ + a){-qr + u, + MVuy - MUvy) = 0. (3.40) 

or 

At X = —L (from j = 4, (3.39)): 

aHq + (^ + a){qr + Ur + MVuy — MUvy) = 0. (3-41) 

OT 

Here the first two terms of the large a expansion of the exact condition are matched 
as is the cr 0 limit. The parameter a could be chosen to optimize the approximation. 
Although we have, in these conditions, been careful to capture the leading order behavior 
as cr oo, we note that terms of the order vajM were earlier neglected. 

So far we have not investigated the well-posedness of the resulting IBVPs, nor carried 
out numerical experiments. In one space dimension, however, we h ave imp lemented 
similar boundary conditions. We note that the troublesome symbol, P, reduces 

to a in one dimension, and so requires no approximation. Therefore, the computations 
only test the accuracy of the small M and small v approximations of the exact conditions. 
Though a complete description of these experiments has appeared elsewhere [7], it is 
worthwhile to show a typical example. Figure 1 contains graphs of q and u for an initial 
pressure pulse. We see the pulse break up into left- and right-moving sound waves which 
propagate through the boundary with no visible distortion. 

In these computations we have used a reasonably fine mesh, 2401 points, and a CFL 
number (based on the sound speed) of .2. The underlying velocity field is in the positive x 
direction, so the left boundary is an inflow boundary and the right an outflow boundary. 
The boundary conditions are then given by: 

u + q = Ux — Ur = 0, X = — 1/2; (3.42) 

u — q = 0, X = 1/2. (3.43) 

These conditions correspond to the specialization of the j = 3,4,5 conditions to the 
one dimensional case. Second order finite differences were used in both the interior and 
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at the boundary. The additional numerical boundary condition sX x — L was second 
order extrapolation of the outgoing characteristic variable ti + g. Full details of these 
computations, including a number of other cases, comparisons with conditions derived 
from energy arguments, as well as higher order (in M) approximations, can be found in 

[7]. 



Figure 1: Plot of q and ix, M = .1, i/ = .01. 


It is interesting to compare these conditions to others which have appeared in the 
literature. Gustafsson and Stoor [5] proposed conditions based on energy arguments. 
Their purpose was to solve slightly compressible flow problems in the absence of sound 
waves, avoiding boundary layers at inflow. We do not expect their conditions to be 
accurate if sound waves are present, and our experiments in one space dimension confirm 
this. More in the spirit of this paper is Halpern’s study of conditions for incompletely 
parabolic perturbations of hyperbolic systems [8]. The philosophy is essentially the same; 
derive expressions for exact conditions and approximate using the smallness of some 
parameters. The small parameters used in that study are u and the tangential wave 
number. (In our case k.) This results in a somewhat different set of boundary relations. 
Our construction has the advantage of allowing k = 0(1), as the assumption (kfa) <C 1 
may be difficult to justify. However, the conditions in [8] do not require M <C 1. 

It is also interesting to compare the derived boundary conditions for the slightly 
compressible model with boundary conditions derived by the same technique for the 
linearized incompressible equations: 

ut + Uu:c + Vuy + px = lyAu, (3.44) 
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Vt + Uvx + Vvy +Py = vAv, 


(3.45) 


Ux + Vy = 0. (3.46) 

Proceeding as above and assuming 0 < u « 1, one obtains the following approximations 
to exact boundary conditions; 

At a: = —L: 

Hp + ut + Uux + Vuy = 0, (3.47) 

Py + Uuy + Ut + Vvy = 0. (3.48) 

At X = L: 

— Hp + ut + Uux + Vuy = 0, (3.49) 

Py + UVx + Vt + VVy = 0. (3.50) 

If we formally set <j = ^P> ^ ~ Mt and use Ux — — Vy then (3.45) becomes (3.38) 
except for the g^-term; (3.46) becomes exactly (3.33); (3.47) becomes (3.37) except for 
the -g^-term; and (3.48) becomes exactly (3.32). There is no boundary condition for the 
incompressible equations corresponding to (3.34). Recall that the inflow condition (3.34) 
comes from the eigenvalue As with large negative real p art (see (3.22)). We note that 
this correspondence requires that our approximation to y/a'^ + approach |A:1 as a — >■ 0. 

To summarize, we’ve derived a nonlinear model system for the study of low Mach 
number flows with sound waves present and systematically derived approximate boundary 
conditions at inflow and outflow for linearizations of the model system about a uniform 
flow. The resulting equations display typical features of accurate conditions for the 
incompressible Navier-Stokes equations combined with a standard radiation condition 
for the wave equation. The latter, however, must be approximated so that long time 
accuracy is obtained. A complete study of the proposed conditions, including numerical 
experiments and analyses of well-posedness and the incompressible limit, are underway 
and will appear elsewhere. 
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